@Suryakanth I tr ied again and noticed that the complaint about an invalid boundary condition is simply due to the fact that tau is not given any value!

So I arbitrarily gave it the value 1.234, as you will see in my worksheet below.

Now what comes up is that the initial Newton integration is not converging.

I would call that some progress.

dust_revised_Solveable.mw

@Suryakanth Sorry to say, but despite several attempts, I'm not getting anywhere.

One attempt was to rewrite the system as follows:

{
F(xi)*diff(diff(F(xi), xi), xi) - diff(F(xi), xi)^2 - 2.200000000*diff(F(xi), xi) - 0.6000000000*xi*diff(diff(F(xi), xi),
 
xi) + 1.000000000*diff(f(xi), xi) = 0, F(xi)*diff(theta_p(xi), xi) - 2*diff(F(xi), xi)*theta_p(xi) - 2.600000000*theta_p(xi) - 0.6000000000*xi*diff(theta_p(xi), xi) + 0.8000000000*theta(xi) = 0,

 (1.058612770 + 1.333333333*(1 + 0.3*theta(xi))^3)*diff(diff(theta(xi), xi), xi) + 1.20*(1 + 0.3*theta(xi))^2*diff(theta(xi), xi)^2 - 14.52406700*theta(xi) - 4.708022333*xi*diff(theta(xi), xi) - 15.69340778*diff(f(xi), xi)*theta(xi) + 7.846703888*f(xi)*diff(theta(xi), xi) + 0.4000000000*theta_p(xi) = 0, 

1.051535713*diff(diff(diff(f(xi), xi), xi), xi) + 1.381281821*f(xi)*diff(diff(f(xi), xi), xi) - 1.381281821*diff(f(xi), xi)^2 - 2.157538185*diff(f(xi), xi) - 0.8287690926*xi*diff(diff(f(xi), xi), xi) - 0.9702972362*M*diff(f(xi), xi) + 0.5000000000*diff(F(xi), xi) = 0
};

I wasn't succesful anyway.

@Suryakanth That works very well!

@Suryakanth What you posted cannot be used in Maple.
You must have done it in Maple originally, certainly.

So upload a worksheet as you uploaded before by using the green up-arrow in the editor.

@Suryakanth  Please take a look at my corrected worksheet below.
One change is that my former local gamma, has been replaced by unprotect(gamma); 

that should work for you too.

Trying to convert OdeSys into a system of first order fails as is repeated in my title above.
dust_revised.mw

@Suryakanth I'm saying you need 9. You say 8.

A comment to my worksheet, which is made in a recent version of Maple, while you have Maple 2018.

I doubt that your Maple version supports beginning like this

restart;
local gamma;

etc......

My guess is that your Maple version complains about assigning to gamma, 

which in Maple stands for Euler's constant (see ? gamma).
Thus you should keep the name gamma_p all over.

I get the message: Too many boundary conditions, need 9, got 12.

See attached worksheet,dust_phase_work_error_BD.mw

I can only see page 1. All the other report the internal server error you mention.
 I use Windows 11.

@nm I don't know the answer to your questions.
I did just now, however, tried odetest on the two returned solutions:

restart;
ode:=diff(y(x),x) = (-y(x)^2+4*a*x)^2/y(x); 
infolevel[dsolve]:=5;
dsolve(ode,y(x), singsol=all);
sol:=%;
odetest(sol[1],ode); #0
odetest(sol[2],ode); #0

It has always been my understanding that dsolve never returns approximate results.

restart;
ode:=diff(y(x),x) = (-y(x)^2+4*a*x)^2/y(x); 
infolevel[dsolve]:=5;
dsolve(ode,y(x), singsol=all);

you will see that Kovacic's algorithm is used succesfully.

simplify(eval((lhs-rhs)(ode),sol)); #0

@dharr The very same is accompliced with subs:

sqrt(1-cos(x)^2);

subs(cos(x)^2=1-sin(x)^2,%);

@nm No what I meant was simply that it gave a correct result if x is real, no more.

## Surely we have for x complex (including x real of course) that:

sqrt(1-cos(x)^2) = sqrt( sin(x)^2 )

## where both square roots are the the principle root.

So this relation holds in any case, but sqrt( sin(x)^2 ) is not the same as abs(sin(x)) for complex x, as e.g. x = 1+ I would show, since clearly abs(z) is real for any z.

This is as far as I get:

A:=sqrt(1-cos(x)^2);
B:=sqrt(sin(x)^2);
simplify(A) assuming x::complex;
%-B;
simplify(%); # 0

This works:

simplify(A) assuming x::real;
# Compare to:
plot(A,x=-Pi..Pi);

This returns NULL:

solve(sol,y) assuming y>1;

and so does this:

solve(sol,y) assuming y<1;